Revealing Majorana Fermion states in a superfluid of cold atoms subject to a 

harmonic potential 
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We here explore Majorana Fermion states in an s-wave superfluid of cold atoms in the presence 
of spin-orbital coupling and an additional harmonic potential. The superfluid boundary is induced 
by a harmonic trap. Two locally separated Majorana Fermion states are revealed numerically based 
on the self-consistent Bogoliubov-de Gennes equations. The local density of states are calculated, 
through which the signatures of Majorana excitations may be indicated experimentally. 
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The Majorana Fermion (MF), a particle that is its own 
antiparticle, was originally proposed by Majorana many 
years ago For the past decade, the MFs have at- 
tracted tremendous attentions due to their exotic prop- 
erties. Especially manipulating MFs is one promising 
way to realize the non-Abelian statistics that may have 
potential applications in topological quantum computa- 
tion 0. For condensed matter systems, the MFs were 
first predicted to emerge in the Pfaffian fractional quan- 
tum Hall state @. Later it was proposed that the MFs 
exist in the p + ip superconducting system associated 
with the zero energy bound states in the vortex core 0- 
6] . While searching for the spin-triplet p-wave supercon- 
ducting materials is of a great challenge. Besides the real 
material, the p + ip superfluid (SF) pairing state was pro- 
posed to be artificially created in the cold atom systems 
with the p-wave Feshbach resonance Q ■ However, p-wave 
molecules are not stable so that the proposal is also diffi- 
cult to implement [8[ . Recently, more attention has been 
paid to topological systems with the spin-orbital interac- 
tion 0-16]. This kind of system with an s-wave pairing 
may be equivalent to the p + ip superconductor 0-[Hl ■ 
Experimentally, tremendous efforts have been made to 
search for MFs. Very recently, several groups have re- 
ported the signatures of MFs in various systems fl7l - l2"il ] , 
but a definite demonstration for the existence of MFs is 
still awaited. 

It has been indicated that both the spin-orbital cou- 
pling and the s-wave pairing SF state can be artificially 
created in cold atom systems. The spin orbital cou- 
pling can be generated through spatially varying laser 
fields (l5l . I22I . [23 | . The s-wave pairing state is much more 
stable than the p-wave one and has been realized [24j . 
Besides the above conditions, the realization of MFs re- 
quire the presence of the Zeeman field and vortex. A 
feasible scheme to realize a controllable Zeeman field in 
cold atoms was proposed in Ref. [15] . In the present work, 
instead of the realization of MFs in the vortex core, we 
elaborate for the first time that the MF states may also 
be induced in a confined harmonic potential, noting that 
the position and motion of the confined potential may be 
controlled more easily in the cold atoms. Based on the 



Bogoliubov-de Gennes (BdG) equations, the SF order pa- 
rameters are calculated self-consistently. Our numerical 
results illustrate the existence of the MF states. The 
distinct features of the MF states are discussed in detail. 

We start from a Hamiltonian including the spin-orbital 
coupling and the SF pairing term, which reads 



H = H t + H A , 



(1) 



Here H t includes the tight-banding part of the model 
with the spin-orbital interaction and the spin polarization 
term, given by [25| 



H t = EM 
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with ipi = (V'if, 4 ! il) T j where a n are the identity (n=0) 
and Pauli matrix (n = 1,2,3), respectively, h represents 
an effective Zeeman field. While H\ is the SF pairing 
term, expressed as 



(3) 



The above Hamiltonian can be diagonalized by solving 
the BdG equations self-consistently, 



H t (r) A(r)a 3 
A*(r)a 3 -a 2 H* t (v)a 2 



*™(r) = E n # n (r) . (4) 



For the N x N lattice size, the dimensionality of the 
Hamiltonian matrix is 4N x <iN. ^"(r) is a 4Af order 
column vector with ^ n (r) = (u^, u™, , v™, , v^) T . 

The order parameters A(r) are determined self- 
consistently, 

A W = J E«T<: + < t <: ) tanh(^|^), (5) 

n 

with V the pairing strength. 

The on-site particle number with spin a is calculated 
after diagonalizing the BdG Hamiltonian, 
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where f(x) is Fermi distribution function. Then we can 
obtain the on-site particle number ri\ = (n^) + (n^) and 
the site-dependent magnetization mi = l/2((n^) — (ni±}), 
respectively. The local density of states (LDOS) can be 
calculated as: 

p-M = V[| < t | 2 6(E n | wfi | 2 S(E n +u)}. (7) 



The delta function S(x) is taken as S — T/[ir(x 2 + T 2 )], 
with the quasiparticle damping T = 0.01. 

We first analyze theoretically the duality between the 
present model with the p-wave superconductors. For the 
uniform case with Ui = 0, the Hamiltonian can be ex- 
pressed in the momentum space, 



H = H, 



Ik 



(A ke ^k,tV'-k,; + /i.c.), 



k 



(8) 



where H t k is the tight-banding term given by 

H tk = ^/i(k)(V' kt Vkt-4^i4)+^[5(k)4t^H+/i.c.] ! 

k k 

with h(k) = (m + cosfc^ + cosfcj,) and ,g(k) = sinfe^ + 
i sin k y . Defining two spinless quasiparticle operators 

Ck = 72(^kT + ei0k ^-H) and d k = jsWlt-e***^), 
the Hamiltonian can be rewritten in the spinless repre- 
sentation: 

H = J2 A k C t (k) ( x C(k) + J2 h(k)CHk)a 3 C(k) 

k k 

+ ^[ 5 (k)e i ^Ct(k)Ct(-k)+/i.c], (10) 

k 

with C(k) = [(d(k), c(k)] T . One can find that the sys- 
tem is dual to the p + ip system with the pairing term 
A — sin h x -\- isinky. Moreover, the spin-orbital inter- 
action is transformed to the pairing term. This dual 
transformation makes it possible to realize the MFs in 
the spin-orbital system. 

The H t in Eq.(2) with t/i = is an effective model to 



describe the topological insulator [25[. The topological 
properties are determined by h with the non-zero Chern 
number for <| h |< 2. Fig. 1(a) illustrates the chiral 
edge states of the model by considering the open bound- 
ary condition along x direction (L x = 50) and periodic 
boundary condition along y-direction with h = —0.2. As 
is seen, the energy gap closes at k y = 7r with the zero 
energy mode being localized at the edge (i x — 0,50), 
suggesting the existence of the edge state. On the other 
hand, if periodic boundary condition is taken for both 
directions, the Hamiltonian can be expressed as 2 x 2 
matrix in the momentum space [Eq.(9)]. The two energy 
bands, obtained by diagonalizing the above 2x2 matrix, 
are plotted in Fig. 1(b), representing the energy spectrum 
in the bulk. There are two bands with the gap 0.4 and 
the energy bandwidth 2 (e 3 [-2.2,-0.2] U [0.2,2.2]). If 
the on-site potential Ui with 0.2 < Ui < 2.2 is chosen, the 




FIG. 1: (a) The energy spectrum of the Hamiltonian Eq. 
(2) with an open boundary condition along the ^-direction, 
(b) The energy band in the momentum space with periodic 
condition. 



Fermi energy would cross the lower energy band and the 
system becomes metal-like. Then, the SF pairing might 
occur when an additional pairing potential term is taken 
into account. 

The presence of the vortex is likely a key ingredient but 
a great challenge for probing the MFs and non-Abelian 
statistics. Usually, the vortex is induced by the magnetic 
field, characterized as a 27r winding of the phase and 
the vanishing of the order parameter in its center. In 
the present work, we will not consider the field-induced 
vortex, rather take into account an harmonic potential 
in the form of Ui = Uq | (rj — r ) | 2 . We can conclude 
from the band structure shown in Fig. 1(b) that the SF 
order parameter is nearly zero in the trap center and 
increases when away from the center. It will reach the 
maximum value when the potential Ui crosses the energy 
band. Two SF boundaries are expected to exist when 
the on-site potential crosses the band edges. The MF 
excitation may appear at the SF boundaries. 

Now let us illustrate numerically the existence of the 
MF states based on the above proposal. In the following, 
we use h = —0.2, Uq = 0.01 and the pairing potential 
V = 2. The numerical calculation is performed on a 36 x 
36 square lattice with the periodic boundary condition. 
The trapping center locates at ro = (18, 18). 

The self-consistent results of the BdG equations are 
presented in Fig. 2. Fig. 2(a) displays the spatial distri- 
bution of the gap magnitudes and their phases (denoted 
by arrows). The spatial distribution of the particle num- 
ber and the site-dependent magnetization are presented 
in Figs. 2(b) and 2(c), respectively. The two-dimensional 
cut of the gap magnitude, the particle number, and the 
site dependent magnetism are plotted in Fig. 2(d). As 
is seen, the gap is nearly zero at the trap center and 
increases when it is away from the center. It reaches 
the maxima at about | ri — ro |= 11. Then the gap 
decreases and becomes nearly zero at the border. The 
SF region forms a ring with the gap magnitudes being 
nearly isotropic around the ring. The phases of the SF 
order parameters change from to 27r around the ring. 
Thus the self-consistent results for the order parameters 
are somewhat similar to the case of the vortex state for 
superconducting materials in the magnetic field. Actu- 
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FIG. 2: (Color online) (a) The numerical results for the mag- 
nitudes and phases (shown by arrows) of the pairing order 
parameters, (b)and (c) are the particle number and mag- 
netism, respectively, (e) The two dimensional cut of the pair- 
ing magnitude, particle number and magnetism along i y = 18, 
respectively. 



ally, here the presence of the phase winding is due to the 
presence of the Zeeman field, which is also an effective 
vertical magnetic field. 

The competition between the SF and ferromagnetism 
is elucidated in Figs.2(b-d), which is a key point for 
the zero mode discussed below. When | rj — ro |< 4.6 
(Ui < 0.2), the on-site energy is smaller than the mini- 
mum excitation energy, and thus the system behaves as 
a ferromagnetic insulator with the particle number per 
site being fixed to near 1.0. The magnetization reaches a 
maximum and the pairing order parameter is negligibly 
small. As 4.6 <| r ; - r |< 15 (0.2 < [/; < 2.2), the 
Fermi energy crosses the energy band. Therefore, the 
lower band becomes unfilled, so that the average on-site 
particle number and the ferromagnetic order decrease. 
As a result, the SF order shows up. We can also see 
that the evolution of ferromagnetism is non-monotonous, 
namely, it reaches the minimum as the pairing order is 
of maximum and recover to a local maxima value as the 
pairing order decreases, indicating the competition be- 
tween SF order and ferromagnetism. As | ri — ro |> 15, 
the on-site energy is larger than the maximum-excited 
energy of the band. Both the particle number and the 
order parameters decrease to zero. 

As presented in Fig. 2, when the on-site energy is small 
or very large, the particle number is fixed to be 2 or 0, 
and thus an insulating gap should exist. The insulat- 
ing gap disappears when the on-site energy crosses the 
energy band. The system becomes metal-like and the 
SF pairing gap shows up. We can expect two effective 
insulator-SF boundaries inside and outside the SF ring. 
The zero mode is expected to exist at the boundaries. We 
diagonalize the Hamiltonian and all of the 5184 eigen- 
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FIG. 3: The eigen-values for the Hamiltonian matrix. Inset: 
The replot of the rectangle part in the main figure. 



values are plotted in Fig. 3. As is seen, the eigen-values 
are nearly continuous but a kink occurs at the zero en- 
ergy. Enlarging the low energy part in the inset of Fig. 3, 
the eigen-values are discontinuous. Intriguingly, our nu- 
merical calculation reveals the existence of a zero energy 
fermionic mode. This zero energy mode is protected by 
an energy gap about 0.027, so that it is robust and can 
hardly be excited by local perturbations. 

As is known, the creation/annihilation operator of the 
captioned zero mode is of self-hermitian namely, 
there may be two MFs associated with this zero mode, 
with the corresponding operators being given by 71 = 
(C + C+J/V5 and 72 = i(C - G)/y/2, where C+ = 
^^{uia^la + Vio-ipia), which can be obtained from the 
zero energy eigen-vector. In this sense, we are able to 
study numerically the MF states by solving the BdG 
equations. The spatial distribution of the two MF states 
Px,2 with the units of p max are presented in Figs. 4(a) and 
4(b), respectively. Since two MFs at the same location 
shall annihilate to an ordinary fermion, only locally sep- 
arated MFs are meaningful. It should be insightful to 
eliminate the overlap part of the two MF states by look- 
ing into the difference of the two MFs (| p\ — P2 |), as 
presented in Fig. 4(c). The two-dimensional cut of p\ — pi 
along i y = 18 is plotted in Fig. 4(d). 

The two separated MF states may be identified from 
the results presented in Fig. 4. They appear at the SF- 
insulator boundaries with maximum probabilities being 
at the horizontal (vertical) directions for 71 (72). One sig- 
nificant feature revealed by Fig. 4 is the non-local behav- 
ior for both states. The MF states are symmetric about 
the trapping center, and thus both of them should have 
two equivalent maxima. Inside the SF ring, there ex- 
ists a certain extent of overlap between these two states; 
while around outsider boundary, there is almost no over- 
lap between the two states. This overlap behavior is a 
distinct feature and would not occur for the usual MFs 
in the well-separated vortex cores. As is known, for the 
usual vortex bounded MF, the wave function decreases 
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FIG. 4: (Color online) The intensity plot of the spatial distri- 
butions of (a) 71, (b) 72. (c) The difference of the probability 
distributions between 71 and 72. (d) The two-dimensional cut 
of pi — p2 along i y = 18. 




FIG. 5: (Color online) (a) The intensity plot of the zero energy 
LDOS. (b) The two dimensional cut of the zero energy LDOS. 
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exponentially away from the vortex core. For the present 
case, the overlap of the MF states seems more interesting 
and may be controlled by the trap potential. 

It is also insightful to look into the LDOS for investi- 
gating the existence and distribution of the above zero 
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In summary, we have proposed and elaborated a new 
scenario to realize the MF states in the fermionic cold 
atom systems. The numerical results have shown that 
the SF region forms a ring with the MF states being 
mainly around the boundaries of the ring. It is rather 
promising that the presence of MF zero modes may be 
detected experimentally through the LDOS spectra. 

This work was supported by the RGC of Hong Kong 
(Nos. HKU7044/08P and HKU7055/09P) and a CRF of 
Hong Kong. 



Rev. Lett. 107, 026405 (2011). 
[12] Jay D. Sau, Roman M. Lutchyn, Sumanta Tewari, and 

S. Das Sarma, Phys. Rev. Lett. 98, 010506 (2007). 
[13] J. Linder et al, Phys. Rev. Lett. 104, 067001 (2010). 
[14] Jay D. Sau, Sumanta Tewari, Roman M. Lutchyn, Tudor 

D. Stanescu, and S. Das Sarma, Phys. Rev. B 82, 214509 

(2010). 

[15] Shi-Liang Zhu, L.-B. Shao, Z. D. Wang, and L.-M. Duan, 

Phys. Rev. Lett. 106, 100404 (2011). 
[16] Xia-Ji Liu, Lei Jiang, Han Pu, and Hui Hu, Phys. Rev. 

A 85, 021603(R) (2012). 

[17] J. R. Williams et al, larXiv: 1202.23231 

[18] L. P. Rokhinson, X. Liu, and J. K. Furdyna. 

larXiv: 1204.42121 

[19] M. T. Deng et al., larXiv:1204.4130l 

[20] A. Das, Y. Ronen, Y. Most, Y. Oreg, M. Heiblum, and 

H. Shtrikman. larXiv: 1205.70731 
[21] V. Mourik et al., Science, 336, 1003 (2012). 
[22] G. Juzeliunas and P. Ohberg, Phys. Rev. Lett. 93, 033602 

(2004); J. Ruseckas et al., ibid. 95, 010404 (2005); P. 



■5 



Zhang et al, Eur. Phys. J. D 36, 229 (2005). 

[23] S. L. Zhu et al, Phys. Rev. Lett. 97, 240401 (2006); 
X.J. Liu et al, ibid. 98, 026602 (2007); T.D. Stanescu, 
B. Anderson, and V. Galitski, Phys. Rev. A 78, 023616 
(2008); C. Zhang, Phys. Rev. A 82, 021607(R) (2010). 

[24] T. Bourdel et al, Phys. Rev. Lett. 93, 050401 (2004); J. 



K. Chin et al, Nature (London) 443, 961 (2006). 
[25] X. L. Qi, Y. S. Wu, and S. C. Zhang, Phys. Rev. B 74, 
045125 (2006); Xiao-Liang Qi, Taylor L. Hughes, and 
Shou-Cheng Zhang, Phys. Rev B 78, 195424 (2008). 



